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ABSTRACT 

Various formulations of smooth-particle hydrodynamics (SPH) have been proposed, intended to re- 
solve certain difficulties in the treatment of fluid mixing instabilities. Most have involved changes to the 
algorithm which either introduce artificial correction terms or violate what is arguably the greatest ad- 
vantage of SPH over other methods: manifest conservation of energy, entropy, momentum, and angular 
momentum. Here, we show how a class of alternative SPH equations of motion (EOM) can be derived 
self-consistently from a discrete particle Lagrangian - guaranteeing manifest conservation - in a man- 
ner which tremendously improves treatment of these instabilities and contact discontinuities. Saitoh & 
Makino recently noted that the volume element used to discretize the EOM does not need to explicitly in- 
voke the mass density (as in the "standard" approach); we show how this insight, and the resulting degree 
of freedom, can be incorporated into the rigorous Lagrangian formulation that retains ideal conservation 
properties and includes the "V/z" terms that account for variable smoothing lengths. We derive a gen- 
eral EOM for any choice of volume element (particle "weights") and method of determining smoothing 
lengths. We then specify this to a "pressure-entropy formulation" which resolves problems in the tradi- 
tional treatment of fluid interfaces. Implementing this in a new version of the GADGET code, we show it 
leads to good performance in mixing experiments (e.g. Kelvin-Helmholtz & "blob" tests). And conserva- 
tion is maintained even in strong shock/blastwave tests, where formulations without manifest conservation 
produce large errors. This also improves the treatment of sub-sonic turbulence, and lessens the need for 
large kernel particle numbers. The code changes are trivial and entail no additional numerical expense. 
This provides a general framework for self-consistent derivation of different "flavors" of SPH. 
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1 INTRODUCTION 

Smoothed particle hydrodynamics (SPH) is a method for solving 
the equations of hydrodynamics (in which Lagrangian discretized 
mass elements are followed; |HIcyl[T977| |Gingold & Monaghan] 
|1977| > which has found widespread application in astrophysical 
simulations and a range of other fields as well (for recent reviews, 
see |Rosswog|2009[ |Spr?ngel|20 10| |Price|20 1 2b) . 

The popularity of SPH owes to a number of properties: com- 
pared to many other methods, it is numerically very robust (sta- 
ble), trivially allows the tracing of individual fluid elements (La- 
grangian), automatically produces improved resolution in high- 
density regions without the need for any ad-hoc pre-specified "re- 
finement" criteria (inherently adaptive), is Galilean-invariant, cou- 
ples properly and conservatively to N-body gravity schemes, ex- 
actly solves the particle continuity equationrl and has excellent 
conservation properties. The latter character stems from the fact 
that - unlike Eulerian grid methods - the SPH equations of motion 
(EOM) can be rigorously and exactly derived from a discretized 
particle Lagrangian, in a manner that guarantees manifest and si- 
multaneous conservation of energy, entropy, linear momentum, and 
angular momentum (Springel & Hernquist 2002, henceforth S02). 

However, there has been considerable discussion in the liter- 
ature regarding the accuracy with which the most common SPH 
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1 This is the continuity equation for a discretized particle field. Exactly 
solving the continuity equation for a continuous fluid, of course, requires in- 
finite resolution or infinite ability to distort the Lagrangian particle "shape." 



algorithms capture certain fluid mixing processes (particularly the 
Kelvin-Helmholtz instability; see e.g. |MorrIsp 996 Dilts] |1999[ 
IRitchie & Thomas | [200T] |Marri & White||2003HOkamoto et~ 
|2003| |Agertz et aT1|2007| >. Comparison between SPH and Eule- 
rian (grid) methods shows that while agreement is quite good for 
super-sonic flows, strong shock problems, and regimes with exter- 
nal forcing (e.g. gravity); "standard" SPH appears to suppress mix- 
ing in sub-sonic, thermal pressure-dominated regimes associated 
with contact discontinuities ( Kitsionas et al. 2009 , Price & Feder- 



rath|2010| [Bauer & Springel|2012||Sijacki et al.|2012H 2 |The rea- 
son is, in part, that in standard SPH the kernel-smoothed density 
enters the EOM, and so behaves incorrectly near contact disconti- 
nuities (introducing an artificial "surface tension"-like term) where 
the density is not differentiable. 

A variety of "flavors" (alternative formulations of the EOM or 
kernel estimators) of SPH have been proposed which remedy this 
(see above and Monaghan 1997; Ritchie & Thomas 2001 



2008 [Wadsley et al. 2008 ; Read et al.|2010|[Read & Hayfield|2012 



— ■ ■■ ■ ' 

Garcfa-Senz et al.|2012p . These approaches share an es- 
sential common principle, namely recognizing that the pressure at 
contact discontinuities must be single-valued (effectively removing 



2 In fairness, we should emphasize that it has long been well-known that 
Eulerian grid codes, on the other hand, err on the side of over-mixing (es- 
pecially when resolution is limited), and in fact this problem actually moti- 
vated some of the SPH work discussed above. This may, however, be reme- 
died in moving-mesh approaches (though further study is needed; see e.g. 
|Springel|20iO) . 
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the surface tension term). Some of these show great promise. How- 
ever, many (though not all) of these formulations either introduce 
additional (potentially unphysical) dissipation terms and/or explic- 
itly violate the manifest conservation and continuity solutions de- 
scribed above - perhaps the greatest advantages of SPH. This can 
lead to severe errors in problems with strong shocks or high-Mach 
number flows, limited resolution, or much larger gradients between 
phase boundaries (J. Read, private communication; see also the dis- 
cussion in [PrTcel|2012b[ |Read & Hayfieldl[20T2l |AbeIp)TT) . All 
of these regimes are inevitable in most astrophysically interesting 
problems. 

Recently however, Saitoh & Makino] {2012) (henceforth 
SM12) pointed out that the essential results of most of these fla- 
vors can be derived self-consistently in a manner that does properly 
conserve energy. The key insight is that the "problematic" inclusion 
of the density in the EOM (as opposed to some continuous prop- 
erty near contact discontinuities) arises because of the ultimately 
arbitrary choice of how to discretize the SPH volume element (typ- 
ically chosen to be ~ m,/p,). Beginning with an alterative choice of 
volume element, one can in fact consistently derive a conservative 
EOM. They propose a specific form of the volume element involv- 
ing internal energy and pressure, and show that this eliminates the 
surface tension term and resolves many problems of mixing near 
contact discontinuities. 

In this paper, we develop this approach to provide a rigorous, 
conservative, Lagrangian basis for the formulation of alternative 
"flavors" of SPH, and show that this can robustly resolve certain 
issues in mixing. Although the EOM derived in SMI 2 conserves 
energy, it was derived from an ad-hoc discretization of the hydro- 
dynamic equations, not the discrete particle Lagrangian. As such 
it cannot guarantee simultaneous conservation of energy and en- 
tropy (as well as momentum and angular momentum). And the 
EOM they derive is conservative only for constant SPH smooth- 
ing lengths (in time and space); to allow for adaptive smoothing 
(another major motivation for SPH), it is necessary to derive the 
"V/i" terms which account for their variations. This links the vol- 
ume elements used for smoothing in a manner that necessitates a 
Lagrangian derivation. And their derivation depends on explicitly 
evolving the particle internal energy; there are a number of advan- 
tages to adopting entropy-based formulations of the SPH equations 
instead. 

We show here that - allowing for a different initial choice 
of which thermodynamic volume variable is discretized - an en- 
tire extensible class of SPH algorithms can be derived from the 
discrete particle Lagrangian, and write a general EOM for these 
methods (Eq. [12] our key result). We derive specific "pressure- 
energy" (Eq. |18[l and "pressure-entropy" (Eq. |21[> formulations of 
the EOM, motivated by the approaches above that endeavor to en- 
force single-valued SPH pressures near contact discontinuities. We 
consider these methods in a wide range of idealized and more com- 
plex test problems, and show that they simultaneously maintain 
manifest conservation while tremendously improving the treatment 
of contact discontinuities and fluid mixing processes. 

2 THE SPH LAGRANGIAN & EQUATIONS OF MOTION 
2.1 A Fully General Derivation 

Following S02, note that the SPH equations can be derived self- 
consistently from the discrete particle Lagrangian 



j N N 

L (q, q) = ~ Yl mi ** ~ S m[ "< 



(i) 



in the independent variables q = (ri, ...,in,h\, namely the 

positions and smoothing lengths of each volume element/particle, 
and the internal energy per unit mass u. If the smoothing lengths 
h are constant, then the only independent variables are the r, and 
the equations of motion follow from d(dL/ 'dqi) / ' (it = dL/dqt. We 
require the derivatives of the recalling that this is at constant 
entropy, so du = — (P/m) dV, we have: 



dm\ _P dAVj 
dqi \ a m, dqi 



(2) 



where p is the pressure and AVj is some estimator of the particle 
"volume." 

We clearly require some thermodynamic variable to determine 
P; we can choose "which" to follow, for example internal energy u 
or entropy A. For a gas which is polytropic under adiabatic evolu- 
tion (we consider more general cases below), if we follow particle- 
carried a,-, then self-consistency with the thermodynamic equation 
above requires that we define the pressure in the above equation 
as P, = (7 — 1 ) Uj (#tj/AVi); if we follow the entropy A, we must 
define P — Aj (m,/AV,) 7 . 

If h is allowed to vary, then we require some relation which 
makes it differentiable in order to make progress. This usually 
amounts to enforcing some condition on the effective "neighbor 
number" or "mass inside a kernel." We stress that this language is 
somewhat misleading: it is not actually the case that there is exactly 
a certain neighbor number or mass inside the kernel, but rather that 
a continuous relation between h and some local volumetric quantity 
is enforced, for example (47r/3)/jj pi = Mkcmci = W;vVngb (so that 
hi oc pj '^ 3 ). Any such constraint (if continuous) is equally valid: 
motivated by the "effective neighbor number" approach we can de- 
fine the constraint equation 



(q) = -W- hi 



1 

M 



-AU=0 



(3) 



where AVi is some continuous estimator of the "particle volume"; 
e.g. for AVi = nn/ pi, we recover the approximate "mass in kernel" 
constraint. 

We stress that AV does not need to be the same as AV; one is 
the effective volume used to evolve the thermodynamics, the other 
is simply any continuous function used to define the hi, so as to 
make them differentiable and thus allow us to include the appropri- 
ate "V/i" terms in the EOM. 

The equations of motion can then be determined from 



d dL dL 

dt dj, ~ % = J 8q~ 



(4) 



where A; are the Lagrange multipliers. The second half of these 
equations (qt — hi) lead to the Lagrange multipliers 



3 Pi AV/ 1 , 

x ' = -^hT^ 



(5) 



hi dAVil, In <9AV,l-' 



SZiU- l^lll I rfil 
3 AVi dh, I 3 AV,- dh, \ K ' 

Inserting this into the first half of the equations gives the EOM 



dv, 

m, — = 
dt 



1\ 

^^[ViAVj + ^jViAVj] 



(7) 



7=1 



Now we require some way of defining "volumes." In SPH 
this is done with respect to the kernel sum: for any particle-carried 
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scalar value x,-, the x-weighted volume average y; = y can be con- 
structed as 



(8) 



7=1 



where WijQii) = WijQr\ij/hi) is the smoothing kernel (discussed 
further below) as a function of \r\tj = |r, — Tj\. The corresponding 
x-weighted volume element can then be defined as AVS = x,/y,. 
Note that for the choice x, = m,, we recover the familiar "standard" 
SPH choices of y, = p\ and AVJ = m;/ pi\ but any other choice of 
Xi is in principle equally valid. For any well-behaved kernel, this 
makes the AVi fully differentiable functions of r and h: 

dAVi _ x, dy, 
dhj yj dhj ' 



V,Ay j = -4V,-y ; - 



(9) 



with 



dyt 
dh 



£f( 3 W + f 



7=1 



dW(\r\/h) I 
0(M/A) La 



Putting all of this together, we obtain 



dv, 
>m — = 
dt 



Y^x,xj [^fijVMjih,) + / V.ll {h :| 



7=1 



Pi 

Xj V3y, art;/ L 3y, dhji 



hi dyr 



(10) 



(11) 



(12) 



(13) 



Because the "potential" (thermal energy) in the Lagrangian 
here depends only on coordinate differences and is rotationally 
symmetric, the pair-wise force in Eq. [T2] is automatically anti- 
symmetric (i.e. obeys Newton's third law); energy, entropy, mo- 
mentum, and angular momentum are all manifestly conserved, pro- 
vided that smoothing lengths are adjusted to ensure the appropri- 
ate (f> constraint (it is straightforward to verify this explicitly). The 
so-called "Vft" terms, which depend on derivatives in time and 
space of the smoothing lengths, are implicitly included to all or- 
ders, via the / terms. The final EOM for any choice of x, and x, 
involves essentially identical information, constraints, and cost. In 
other words, because all of the formulations we consider are sim- 
ply replacing the choice of particle-carried scalar in Eq. [12] they 
involve identical computational expense for otherwise equal gas 
conditions. 

2.2 Formulations of the SPH Equations 

2.2.1 Density-Entropy Formulation 

If we take x; = x, = m, (giving y, = y t — p,, the kernel-averaged 
mass density estimator, and volume estimator AV; = mt/ pi), and 
follow the entropy A,- (giving P, = Pj = AipJ), we obtain the EOM 
from S02: 



dv; 
df 



7=1 



Pj 



[fiAipy-^iWijihi) +fjAjp]- 2 V i W ij (hj)] 



7=1 



*-[ 



Note that for adiabatic evolution, we require no energy equa- 
tion, since entropy is followed; for a specific energy defined as 
u = Pi/[(j— I) pi] = (7 — 1) _I Aip]~ l , energy conservation is 
manifest from the EOM above. 

As discussed in § [TJ this formulation is known to have 
trouble treating certain contact discontinuities. Because the only 
volumetric quantity that enters is the density, this fails when 
the densities are no longer differentiable, even when pressure 
is smooth/constant. Specifically, consider a contact discontinuity, 
pi c\ = p2c\ (where quantities T and '2' are on either side of 
the discontinuity). As we approach the discontinuity, the kernel- 
estimated density must trend to some average because it spherically 
averages over both "sides," p — > (p), but the particle-carried sound 
speeds c remain distinct, so the pressure is now multi-valued, with 
a 'pressure blip' of magnitude ~ (p m ax/pmin)firue appearing. This 
has a gradient across the central smoothing length, causing an ar- 
tificial, repulsive "surface tension" force that suppresses interpene- 
tration across the discontinuity. 

2.2.2 Pressure-Energy Formulation 

Instead consider Xt = Xt = (7 — 1 ) I/; = (7 — 1) rm u, , proportional 
to the particle internal energy. Now, y, = y"; is by definition a di- 
rect kernel-averaged pressure estimator y, = Pi, and AV; = (7 — 
i)mjUj/Pj. This means that the pressure itself is now the directly 
kernel-averaged quantity entering the EOM and is therefore al- 
ways single- valued. So long as the pressure is smooth/differentiable 
(regardless of how the density varies), the EOM should be well- 
behaved. For this choice, we obtain)^] 

^ = -£(7- ^f m i »■ "7 [j. VWiihi) + & VtWijihj)] 

j=l ' J 

- N 

f,= [ l + JpMX^ fl = £(7-l)«^WX*») (15) 
' ' 7=1 

Recall, we now need to evolve the energy explicitly. From 
Eq.[2] du/dt — — (P/m) (dAV /df), and the evolution of the volume 
element AV just follows from the Lagrangian continuity equation, 
so we obtain 



~d7 



= £(7-1 



7=1 



2 fi 

) rrij UiUj -= (vj - \j) ■ VjWjj(hi) 

Pi 



(16) 



It is straightforward to see that this guarantees explicit energy con- 
servation; entropy conservation is implicit and also easily verified. 

There are, however, significant drawbacks to the choice of 
Xi = Xi (AV — AV), in which case we are implicitly defining h such 
that h] oc iij/Pj. This is, by itself, perfectly valid and is easily solved 
by the same bisector method as in the "standard" (AV = m/p) for- 
mulation. However, in practice, the particle values vary much 
more widely than the m,. This leads to some potential problems. 
First, if there is large variation in u,, convergence in /j, can become 



3 Eq.^Jis similar to the EOM derived in SM12, itself identical to that de- 
rived earlier in Ritchie & Thomas 12001 1 from purely heuristic arguments. 
These do, after all, motivate the derivation here. However there are two key 
differences. First, we derive and include the Vft terms (/, = 1 in SM12), 
necessary for conservation if h varies. Second, the \7W terms enter dif- 
ferently (with different multipliers and indices). This stems from the La- 
grangian derivation and is necessary - even for constant h - for the EOM 
to simultaneously conserve energy and entropy (i.e. to properly advect the 
thermodynamic volume element). 
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quite expensive. Second, again under circumstances with large dis- 
order, the required hi can become very large, leading to an effective 
loss of resolution. Third and most problematic, under some circum- 
stances the constraint cj> can have multiple solutions; in this case if 
hi "jumps," it is no longer continuously differentiable, and so exact 
energy conservation is broken. 

An obvious alternative is to use x = 1, i.e. AVi = 1 /«,, where 
hi is the "particle number density" 



m = yf(x = 1) 



7=1 



(17) 



This restores the effective "number of neighbors" criterion for hi, 
and is always well-behaved since all particles are weighted equally. 
If we do this, the EOM become: 



dVi 
dt 



5^(7- ifmjUiUj \&ViWij(hi) + ^ViWij(hj)\ 
OP, 



fa- 



I . u [\ — mw+^r as) 

\3(7— IjUjmjUj oh,/ L 3n, o/j,J 

Note that this is just the previous equation with fj — s> fy\ in other 
words, the EOM are identical up to the "VA" corrections, which 
is what we expect, since the only function of the AV term is to 
determine how the hi evolve. Trivially, then, the energy equation is 
also the same as above but with f, — > fij. 

As discussed in SMI 2, because the volumetric quantity used 
in the EOM here is now directly the kernel-estimated pressure (in- 
stead of the density), this formulation automatically guarantees that 
pressure is single-valued at contact discontinuities, and so removes 
the pressure "blip" and surface tension force. The equations will 
now be well-behaved so long as pressure is smooth. This is true 
by definition in contact discontinuities; it is of course not true at 
shocks, but neither (typically) is the density constant there - so we 
do not lose any desirable behaviors of the density-entropy formu- 
lation. In either case, we require some artificial viscosity to treat 
shocks. 

2.2.3 Pressure-Entropy Formulation 

If we wish to retain a direct kernel-estimate of the pressure entering 
the EOM, but formulate this in terms of entropy, we must instead 
consider x, = nijA 1 / 1 . In this case|^]we obtain by definition (from 
the consistency requirement for Pi) 



P i =y]=[j2>n J A l /-<W, J (h,)}' 

7=1 



(19) 



If we also assume x, = x,-, the EOM become 



dv, 
dt 



./=! 



mj (AiAj) * [|£ ViWijih) + £gL V,W t j{hjj\ 



fjPj 



h, dP] h 



3=1/7 dh 



-V 

% J 



(20) 



4 This choice of Xi may seem a bit strange, but in fact this is the only self- 
consistent "entropy formulation" which directly evaluates the pressure. If 
we simply substituted uj = pj 1 / (7 — 1) in x, = (7— 1) m, h,, we would 
re-introduce the density p (which we are trying to avoid in this formulation 
of the equations); we could instead define Uj = Aj {Pj/Aj)~<~~ [ / (7 — 1), 
but this involves Pj in its own definition and would require a prohibitively 
expensive iterative solution over all particles every timestep. 



Note however that using Xi = x, implies AVi = m, (Ai/Pi) 1 ^; 
as in the previous section, this can introduce problems of conver- 
gence and diffusion. Therefore, instead consider as before x; = 1 
(AVj = !/«;)• This gives: 



dv, 
dt 



with 



- ]T mj (A,A^ [ ViWij(hi) + f -0i ViWij(hj)] 



fij 



1 



dPi 



l/ 7 



3 A 1 / 7 mj nj 



dh, 



[i+A^r 

V 3h, dhil 



(21) 



(22) 



As in the density-entropy formulation, we explicitly evolve the 
entropy so for adiabatic evolution require no additional evolution 
equation. 

This formulation is very similar to the pressure-energy for- 
mulation, (and has the identical advantages of good behavior at 
contact discontinuities). The only difference is the free choice of 
thermodynamic variable. This formulation trivially conserves en- 
tropy, and manifestly conserves energy to machine differencing ac- 
curacy if constant timesteps are used (the choice of pressure-energy 
or pressure-entropy formulation can lead to some differences when 
adaptive timesteps are used, but we show these are generally small). 
It is largely a matter of convenience and minor computational ex- 
pense which method is preferred. 

When the pressure is smooth and there is good particle order, 
the fij ~ 1 here, which means our choice of how to regularize h is 
unimportant, and no spurious "surface tension" force is introduced. 
For the choice x = 1, the correction terms remain well-behaved 
even if there is large particle disorder in A,-, critical to stability 
in simulations when heating/cooling are included and entropy is 
no longer conserved. Another useful feature here is the following: 
imagine the case where there is large particle disorder so Pj 2> Pi 
and Aj 3> A,. Since the A terms enter as multiplicative pre-factors, 
their difference does not introduce errors into the sum; gradient er- 
rors will arise from differencing the P, terms, but for 7 = 5/3, these 
enter only as P t '^ 5 , so differencing errors are greatly suppressed. 

2.2.4 More General Cases 



In § |2.2.1|2.23| we simplify by assuming the gas obeys a poly- 
tropic equation of state under differential adiabatic compression or 
expansion. We emphasize that this does not exclude the gas under- 
going shocks (in which the entropy and energy change according to 
artificial viscosity), cooling, and/or chemical evolution (additional 
operations to duj in Eq. |2j; these are just handled in an additional, 
separate step or loop each timestep (see an example in § |4.8[ >. 

However, some situations call for more complicated equations 
of state. Consider the case where the pressure of a given particle 
Pi is an arbitrarily complicated (but single-valued) function g of 
the thermodynamic volume element AVi and the local (particle- 
carried) state variables a, = (a,. 1 , a,, 2, ...,«,.„,), so P, = g(AVi, a,). 
The a might include m, and «, or A, , as in our previous examples, 
but also information about the chemical state, radiation field, po- 
sition or velocity, phase, etc, of the gas. Our general form of the 
EOS in Eq. [T2] made no assumption about the equation of state, 
and still holds. The question is how to determine the appropriate 
x, and Xi for a "pressure formulation." This requires any x, such 
that there is a one-to-one mapping between the smoothing kernel 
sum and the pressure (so that V(AV) vanishes when VP does). 
We can ensure this by choosing x, to be the solution to the equation 
g(AV; = Xi, a,) = 1 (i.e. if we were to replace AV, by x,, which 
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we recall typically has different units, we would obtain a dimen- 
sionless P, = 1). We then define AV; = Xi/yt = Xi/^2xjWij(hi) as 
usual, and take P, = g(AVj = Xi/yi, a,) in Eq. 12 We still have full 
freedom to define xf, based on the formulations above, we suggest 
that the simple choice x, = 1 is often most stable. 

It is straightforward to verify that our choices in both the 
pressure-energy (p = (7 — l)i(,m,/AV;) and pressure-entropy 
(P = Aj (nii I AViY') formulations satisfy the condition g(AVi = 
Xt, a,) = 1. These are just special cases of the above with a, = 
(nij, in) or a,- = (m,, A,) (and different g), respectively. 

A more complicated example might be a case with a chemical 
potential, such that m, diii\A = — Pi AAVj + ^ M* d/Vi- where the k are 
different species. If the Nt are constant over a differential adiabatic 
volume change, then this can be treated by any of the above formu- 
lations with the chemistry evolved (under the effects of radiation, 
for example) in a separate step. If the 2V* are functions of the volume 
element AV, however, then we simply have rmdui\A — ¥ — PdAV; 
where P, = P - J2 ^ (dN k /dAVj). So we can use the EOM in 



Eq. 12 with P, — ¥ P, and use the approach above to determine the 
appropriate x, (wifhx, = 1). 

3 ADDITIONAL SIMULATION INGREDIENTS 

3.1 Entropy & Artificial Viscosity 

As is standard in SPH, the algorithm is inherently inviscid and some 
artificial viscosity must be included to properly capture shocks. 
However, we include a more sophisticated treatment of the artificial 
viscosity term (as compared to S02 and "standard" GADGET; which 
follow Gingold & Monaghan 1983 1, following Morris & Monaghan 
l |1997) with a|Balsara ( 1989} switch. This includes a particle-by- 
particle artificial viscosity that grows rapidly in strong shocks and 
rapidly decays away from shocks (to a minimum a = 0.05), reduc- 
ing numerical dissipation by more than an order of magnitude away 
from shocks compared to the previous constant artificial viscosity 
prescription. For detailed comparison of the viscosity algorithms, 
we refer to |Cullen & Dehnen| ( |2TJlO) , 

3.2 Thermodynamic Evolution & Timestep Criteria 

For all problems discussed here, we employ the GADGET adap- 
tive timestep algorithm, which dramatically reduces the computa- 
tional expense for almost all interesting problems (relative to using 



a constant simulation timestep). However, as pointed out in [Saitoh 



& Makino (20091 and developed further in |Durier & Dalla Vec 



chia (20121, in problems with very high Mach number shocks/bulk 



flows, the standard adaptive timestepping can lead to problems if 
particles with long timesteps interact suddenly mid-timestep with 
material evolving on much shorter timesteps. Fortunately this is 
easily remedied, and we do so by implementing a timestep limiter 
identical to that in |Durier~& Dalla Vecchia (2012). At all times, 
any active particle informs its neighbors of its timesteps and none 
are allowed to have a timestep > 4 times that of a neighbor; and 
whenever a timestep is shortened (or energy is injected in feed- 
back) particles are "activated" and forced to return to the timestep 
calculation as soon as possible. 

3.3 Smoothing Kernel 

Our derivation of the EOM allows for an arbitrary choice of SPH 
smoothing kernel W, so long as it is differentiable. This choice can 
have a significant effect on some test problems via its effect on 
the resolution pressure gradient errors (see e.g. Morris 1996j |Dilts| 
1 1999} [Read et al. 2010, and discussion below). We have experi- 
mented with a wide range of possible kernel shapes, following the 



~ 10 discussed in |Fulk & Qumn| ( [T996) and|Hongbin & Xin| (2005 ), 
the triangular kernels proposed in |Read et al.| ( |2010} , and the vari- 
ant Wendland kernels proposed in Dehnen & Aly (2012). However 
our intention here is not to study SPH kernels; for simplicity we 
therefore adopt a standard quintic spline kernel with Nngb = 128 
neighbors in all tests shown here (unless otherwise noted). This is 
the 'optimal' spline kernel suggested in both Hongbin & Xin ( 2005 ) 
and Dehnen & Aly (2012) and has effective resolution equal to a 
cubic spline with 34 neighbors (but significantly higher accuracy). 
In most cases, we obtain similar results using a lower-order cubic 
splin^jwith Nngb = 32; but we discuss where this is not the case. 
We do not see qualitative improvement in the specific tests here us- 
ing yet higher-order kernels and/or increasing neighbor number as 
high as /Vngb ~ 500. 

3.4 Density Estimation in Density-Independent SPH 

An estimate of the density is often required for other calculations 
(e.g. cooling), even if it is not needed for the EOM. This can still 
be directly estimated from the standard SPH kernel, as p, above 
(i.e. a kernel sum with x, = mi). But in principle a density could 
also be inferred or estimated as p,h fa m;/AV; for the thermody- 
namic AV "pressure formulations" (i.e. from the combination of 
the variables P; and A, or ui). However, in "mixed" regions near 
contact discontinuities, this will lead to multi-valued densities in 
neighboring particles (since entropies are still conserved at the par- 
ticle level). This relates to the fact that the two densities represent 
physically distinct quantities. The former, direct kernel pi (x, = mi) 
is simply the volume-average mean mass density at the particle lo- 
cation. The latter (thermodynamically-inferred p r ;,) is the energy 
or entropy-weighted mean mass density (in the pressure-energy 
or pressure-entropy formulations, respectively), averaged over all 
neighbors within the kernel. Because the cooling is typically calcu- 
lated on a particle-by-particle basis, we find the former estimator is 
more stable and appropriate in most applications. However, there 
may be situations where alternative weightings for the density es- 
timator are optimal, and a more complete treatment of mixing may 
require a mechanism for equilibrating entropies inside the kernel 
and generating mixing entropy. 

4 TEST PROBLEMS 

4.1 Strong Sedov-Taylor Blastwaves 

Here we consider an extreme Sedov-Taylor blastwave, with very 
large mach number, designed to be a powerful test of conservation. 
A box of side-length 6kpc is initially filled with 128 3 equal-mass 
particles at constant density n — 0.5 cm~ 3 and temperature 10 K; 
6.78 x 10 46 Jof energy is added to the central 64 particles in a top- 
hat distribution. This triggers a blastwave with initial Mach number 
~ 1000, which we compare at 20Myr, where the shock front should 
be at r ~ 1.19kpc. In this test and below, unless otherwise specified, 
we assume a 7 = 5/3 gas equation of state. 

First we compare the "fully conservative" algorithms we 
derive above. In "standard" SPH (density-entropy formulation; 
Eq,|14}, the analytic solution is reproduced very well up to the SPH 
smoothing/resolution limit. The narrow shock jump is not perfectly 



5 It is well-known that the "standard" Schoen berg| jl946) spline kernels 
become unstable for arbitrarily high neighbor number; the "correct" way 

to increase /Vngb to better sample the kernel and reduce gradient errors is 

1/3 

to move to higher-order kernels of 0(n) fs — 1 + (1.0 — 1.2)/V N q B . Done 
properly, this allows increasing /Vngb an d accuracy without decreasing the 
effective spatial resolution (albeit at additional computational expense). 
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Figure 1. Strong Sedov-Taylor blastwave test. We compare the solution at fixed time and resolution, with different SPH algorithms discussed in the text. We 
plot radial velocity, density, and temperature (for clarity, we plot the median and 1% — 99% interval of particle values binned in radial intervals Ar = 0.01), 
with the analytic solution. We measure the accuracy of energy conservation SE/E = \E(t) — E(t = 0)\/E(t = 0). Left: Standard SPH (density-entropy; 
Xj = Xj = m,); Eq. |14| Center: Lagrangian Pressure-Energy formulation (Xj = (7 — 1) £/;) with particle-number density based smoothing lengths (,v, = 1); 
Eq.^j Right: Lagrangian Pressure-Entropy formulation (jc, = m/A*^ 7 ) with particle-number density based smoothing lengths (x, = 1); Eq. |2l| Conservation 
and accuracy is excellent up to the resolution limits (lack of particles dominates the noise at small radii), in all three algorithms. The pressure formulations 
slightly increase the particle scatter in temperature, but reduce the post-shock "ringing" of the velocity solution. Conservation errors in all three cases are 
overwhelmingly dominated by the adaptive timesteps, not the EOM. 
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Figure 2. As Fig.^ but for SPH equations that are not explicitly conservative. Top: Pressure-Energy formulation without the Vft terms (ft = 1 in Eq. |15) . 
Variations in h now cause order-unity energy conservation errors, leading to the shock being in the wrong position. Middle: Pressure-Entropy formulation 
without the S7h terms. Again order-unity conservation errors appear and the shock evolves incorrectly. Bottom: Non-Lagrangian SPH as considered in |Morris] 
(1996); Abel 12011 1; this algorithm minimizes linear pressure errors and removes the surface tension term, but violates conservation of energy and momentum. 
Conservation errors grow exponentially and dominate the solution. 



resolved, for example, but averaged over the same smoothing, the 
agreement is excellent between analytic solution and mean parti- 
cle values. There is also some scatter in particle properties (most 
notably velocity), but it is generally narrow. The largest deviation 
from the analytic solution is in the post-shock "ringing" in the ve- 
locity field, a well-known effect (which is sensitive to the artifi- 
cial viscosity prescription). Note that the behavior at small radii 
< 0.5 kpc is noisy simply because the extremely low density means 
there are no particles here. As expected, the conservation properties 
are excellent as well; energy is conserved to a part in ~ 10 4 ; simi- 
lar conservation obtains for entropy, energy, momentum, and angu- 
lar momentum. In fact the conservation errors are overwhelmingly 
dominated by the adaptive time-stepping scheme (which violates 
conservation by not always evolving mutually interacting particles 
at the same timesteps), not the formulation of the EOM. If we sim- 



ply force constant, global timesteps, the conservation errors are re- 
duced to machine accuracy (however this comes at great numerical 
expense)]^] 

To compare, the pressure-energy formulation (jr, = (7 — 1) Ut 
as Eq. [T8] using the neighbor number density volume element 
Xj = 1 to define h) agrees very well, and also gives excellent conser- 
vation. The particle noise/scatter is slightly larger, most noticeably 
in the post-shock temperature. This occurs because of some "mix- 
ing" of thermal properties in the kernel average from higher-energy 



6 We do confirm the point discussed in § |3.2| that motivates our timestep- 
ping algorithm: without care in "signaling" when adaptive timesteps are 
taken, conservation errors can be severe. 
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Figure 3. As Fig. [T] with explicitly conservative equations, but different 
choices of .v, (how to regularize the hi). Top: Lagrangian Pressure-Energy 
formulation with Xj = .v, (Eq. |15) . Although conservation is maintained, this 
choice of Jc, leads to enormous numerical diffusion and particle disorder 
("spreading" the shock location over a large radius and affecting the pre- 
shock gas). Bottom: Lagrangian Pressure-Entropy formulation with Jf, = x, 
(Eq.pSJ. Again, this leads to severe diffusion. 



particles. However the post-shock ringing in the velocity solution 
(although still present) is reduced. 

The pressure-entropy formulation (x, = nijA 1 / 1 with x, = 1, 
Eq. \2\\ is essentially identical to the pressure-energy formulation 
here (the slightly better conservation owes to how the sound speeds 
enter the timestep signaling scheme). This is not surprising - the 
two represent essentially identical information but just choose to 
explicitly follow different thermodynamic variables. 

As discussed in § |2.1| the EOM for the pressure-energy and 
pressure-entropy formulations take a simpler form if we assume 
Xj = xt, as in standard SPH, but this can lead to problems. We 
show that here, by considering these implementations (Eq.|15|&|20| 
for the pressure-energy and pressure-entropy formulations, respec- 
tively). Recall, these are still fully conservative Lagrangian formu- 
lations; as a result we still see good energy conservation. However, 
forcing the smoothing lengths to evolve not just with particle num- 
ber density but with local thermodynamic quantities leads to oc- 
casional enormous "super-smoothing" that is both computationally 
expensive and introduces enormous numerical diffusion. We see 
that here, where the upper envelope of particles and even the mean 
are biased in the pre-shock medium (information from the post- 
shock region has clearly been over-smoothed into the regions here). 
Although some mean and post-shock quantities are well-behaved, 
this algorithm has introduced far too much numerical mixing. 

It is also instructive to see what happens if we drop the V/i 
terms in these formations (returning to x, ■ — 1). The EOM will now 
violate manifest energy conservation to the level that the smooth- 
ing lengths vary over an individual smoothing length (in a smooth 
medium, this correction should vanish at infinite resolution, but it 
will never vanish in shocks if variable h are allowed). Indeed we 
now see order-unity energy errors. Although the qualitative solution 
appears similar, the pressure-energy and pressure-entropy solutions 
gain and lose energy, respectively, and so the shock is simply in the 
wrong position. This occurs in other, similar non-conservative for- 
mulations as well, for example see the discussion of problems with 
|Ritchie & Thomas|p00l) & 'OSPH' methods in |Read & Hayfield| 
( |2012| Appendix A). 

Finally, we consider an algorithm which is manifestly non- 
conservative. Specifically, we consider the Morris ( 1996) formula- 
tion of the pressure derivative, where the kernel sum for the EOM is 
over (Pj — Pi)/(pi Pi) VWij. This is very similar (although not iden- 
tical) to the EOM proposed in Abel ( 201 1) as well. It is possible 
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Figure 4. Sod shock tube in three dimensions at time t = 0.1. We show 
the median particle density and pressure along the long x-axis position 
(binned as Fig. [TJ, compared to the reference solution from a piecewise 
parabolic method (PPM) solver. Inset shows the pressure profile binned 
in much smaller intervals around the location of the contact discontinuity 
(x m 0.44), where the pressure should be constant. Top: Standard (density- 
entropy) SPH. The exact solution is generally reproduced well up to SPH 
smoothing, as Fig.[T] but a "pressure blip" with an artificial gradient appears 
near the contact discontinuity. Bottom: Pressure-entropy SPH (Eq.[2T] with 
xi = 1 for the reasons in Fig. [3j. The pressure blip is reduced to particle 
noise-level, while the agreement elsewhere with the PPM result remains. 



to show that this formulation actually eliminates the leading-order 
gradient errors associated with the "standard" SPH EOM (see e.g. 
|Price||20"l2b| >. However, the cost is manifest resolution-level vio- 
lation of conservation (of energy, momentum, and entropy). If we 
adopt our standard initial conditions with this algorithm, we see 
that the conservation errors grow exponentially and quickly swamp 
the real solution. The problem, as described in |Price| ( |2012b} is in 
the non-linear terms that maintain particle order in sph[] 



7 |Abel|(201l) do show that the behavior of Sedov blastwaves in their for- 
mulation is tremendously improved if the initial conditions are modified 
so that the blastwave does not start from a point injection or top-hat parti- 
cle distribution, but from an already-developed smaller blastwave with pre- 
initialized, resolved pressure gradients. However that is not the particular 
test here, and - as they caution - is not often the case in astrophysical sys- 
tems where e.g. early-stage SNe explosions are unresolved. 
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4.2 Sod Shock Tube 



t = 



Fig.[4]shows the results of a standard three-dimensional Sod shock 
tube test. We initialize a period domain with lengths 2, 1/8, 1/8 in 
the x, y, z directions, 5 x 10 4 particles, 7 = 5/3, zero initial veloci- 
ties, densities p = 1 , 0.25 and pressures P = 1 , 0.22 in the left/right 
halves of the x domain. In Fig. [4] we compare the results from 
both standard density-entropy SPH and the pressure-entropy for- 
mulation (with Xi = 1) at time t = 0.1. The median particle values 
agree well with the exact solution in both cases. As expected and 
seen in the Sedov case, the density and pressure discontinuities are 
smoothed by the SPH smoothing and artificial viscosity but oth- 
erwise well-behaved. The pressure discontinuity is slightly more 
smoothed by direct kernel averaging in the pressure-entropy case, 
but the difference is small. 

If we zoom in on the pressure profile near the contact discon- 
tinuity at x ~ 0.44 (where the exact solution is P — constant), we 
see the "pressure blip" discussed in § |2.2.1| appear in the density- 
entropy case. The presence of the smoothed density in the EOM 
leads to an artificial pressure gradient; this occurs over a couple 
of SPH smoothing lengths with a fractional amplitude of ~ 10%. 
However, it is clearly above the particle-noise threshold, and we 
show below that it has significant effects. In the pressure-entropy 
formulation, this is almost completely eliminated (at least reduced 
to the noise level), whether we choose x, = 1 or x, = x ; . 

We have also repeated the ID shock tube tests in SM12; our 
results are generally indistinguishable from their Figs. 2-4. 

4.3 Hydrostatic Equilibrium/Surface Tension Test 



Here we consider a simple test following Cha et al. (2 010| l; |HeB &| 
Springel (2010) and SM12, that allows us to see the consequences 
of the "surface tension" effect discussed in § Q] We initialize a 
two-dimensional fluid in a square of length L = 1 (with period 
boundaries) and constant pressure P = 3.75, polytropic 7 = 5/3, 
and density p — 4 po within a central square of length L = 1/2 and 
p = po = 7/4 outside. We use 256" total particles, though the re- 
sults are similar for as few as ~ 50 2 . 

Fig.[5]shows the resulting system at t — 0, and evolved to t = 3, 
in the "standard" density-entropy formalism (Eq.|14[> and pressure 
entropy-formulation (Eq. |21| it makes little difference for this test 
whether we adopt x, = x, or x, ■ = 1 ). This should be a stable configu- 
ration. But the density-entropy case behaves as a system with phys- 
ical surface tension; a repulsive force appears on either side of the 
contact discontinuity from the "pressure blip," opening the gap in 
the plot which then deforms the square to minimize the surface area 
of the contact discontinuity. It converges to a stable circle after a 
few relaxation oscillations. In the pressure-entropy case, the square 
remains stable for the duration of the runs we consider (t = 50). 
There is a slight "rounding" of the comers, but this occurs quickly 
(? < 1) then stabilizes; it appears to be a direct consequence of the 
smoothing kernel. These results are identical to those in SMI 2, for 
the pressure-energy formulation (with no Vh terms). 

4.4 Kelvin-Helmholtz Instabilities 

We next consider a (three-dimensional) Kelvin-Helmholtz (KH) 
test. The initial conditions are taken from the Wengen multiphas e 



Read et al 



i 2010 i 



test suit^jand described in Agertz et al. 1 2007 
Briefly, in a periodic box with size 256, 256, 16kpc in the x, y, z di- 
rections (centered on 0, 0, 0), respectively, ~ 10 6 equal-mass par- 
ticles are initialized in a cubic lattice, with density, temperature, 




Pressure-Entropy 




Available at http : / /www . ast rosim . net /code/doku .php 



Figure 5. Density (red = 4x blue) in a hydrostatic equilibrium test, with 
uniform pressure and no external forces. Top: Initial condition. Middle: Sys- 
tem evolved to time t = 3 in the pressure-entropy formulation (Eq. |21) . The 
square evolves stably (the small corner rounding stems from the smoothing 
kernel). Bottom: Time t = 3 in standard (density-entropy) SPH. The "pres- 
sure blip" around the contact discontinuity (Fig. [4) manifests as an effective 
surface tension, opening a smoothing-length gap between the two fluids, 
and gradually deforming the square into a circle. 



and x-velocity = p\ , T\ , vi for \y\ < 4 and = p%T%, V2 for \y\ > 4, 
with p2 = 0.5 pi, T2 = 2.0 Ti, V2 = — Vi = 40kms~'. The values 
for T\ are chosen so the sound speed c s ,2 ~ 8 1 V2 1 ; the system has 
constant initial pressure. To trigger instabilities, a sinusoidal veloc- 
ity perturbation is applied to v y near the boundary, with amplitude 
Sv y = 4kms -1 and wavelength A = 128kpc. 

Fig. [6] compares the resulting behavior in the "standard" 
density-entropy formulation of SPH (Eq. \14\ , and the pressure- 
entropy formulation (Eq,|21[> proposed here. Because of the severe 
diffusion seen in the Sedov test associated with the choice of x, = x, 
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Figure 6. Specific entropy map of Kelvin-Helmholtz instabilities at time 
2 times the linear K-H growth timescale tkh (top), and at 8tkh (bottom). 
Left: Standard (density-entropy) SPH. Note the "surface layer" surrounding 
the curls and the breakup into oil-and-water style blobs at later times. Right: 
Pressure-entropy formulation (.?,- = 1). 



for the pressure-entropy formulation, we focus on the formulation 
with the "neighbor number density" volume element used to define 
h, i.e. xt = 1. However, for the sub-sonic, pressure-equilibrium sys- 
tems simulated in this test, the choice of x, makes little difference 
(highlighting the importance of our strong shock tests). We also do 
not further separately consider the pressure-energy formulation, as 
it gives very similar results in our subsequent tests. As expected 
from the discussion in §[T] the density-entropy formulation poorly 
captures the KH instability. The surface tension term forms a sharp 
boundary layer between the two phases, which suppresses all the 
small-wavelength mixing and leads at late times to the breakup of 
the "rolls." This ultimately produces blobs that resemble an oil-and- 
water morphology (a system with physical surface tension). 

In the pressure-entropy formulation, however, the behavior 
is dramatically improved. The long-wavelength mode grows on 
the correct (linear) KH timescale, but we can also plainly see the 
growth of modes at all wavelengths (along the "edge" of the den- 
sity surface) down to the kernel scale. At later times, we resolve 
~ 3 — 5 full "wraps" of the rolls while the standard SPH solution 
breaks up. Most important, there is obviously actual mixing occur- 
ring "inside" the rolls. Direct comparison to high-resolution results 
from grid codes (both fixed-grid and adaptive mesh solutions) as 
well as Godunov SPH and moving-mesh methods, using identical 
initial conditions, show that the result here is quite similar (compare 
e.g. |Readetal.|2010| Fig. 6, or |Murante et al.|2011| Figs. 5-9). 

In Fig. [7] we repeat this experiment with the pressure-entropy 
formulation, but replace our standard treatment of artificial vis- 
cosity (see § |3.1| > with a more simplified treatment (using a con- 
stant artificial viscosity for all particles and times); the latter is 
known to produce significantly more numerical dissipation away 
from shocks. Although it is well-known that this can produce sig- 
nificant differences in some situations (e.g. sub-sonic turbulence 
and rotating shear flows; see |Cullen & Dehnen|2010|[Price & Fed-| 




Figure 7. Comparison of the Kelvin-Helmholtz behavior at t = 12tkh on 
artificial viscosity in the Pressure-Entropy formulation. Left: Our standard, 
by-particle time-dependent artificial viscosity following Morris & Mon- 
|aghan| [l997 Right: The original GADGET simplified constant artificial 
viscosity prescription (as Gingold & Monaghan 1983 with a Balsa ra|1989| 
switch). Although very important for some behaviors, it has little effect on 
KH instabilities within this formulation. 
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Figure 8. Comparison of the Kelvin-Helmholtz behavior at t = 3 tkh on 
the smoothing kernel in the Pressure-Entropy formulation. Top Left: Quar- 
tic core-triangle kernel with neighbor number Nngb = 442 from |Read| 
|et al.|j2010) . Top Right: Cubic core-triangle with N NaB = 128. Bottom 
Left: Quintic spline with Wngb = 128. Bottom Right: Cubic spline with 
%GB = 32. 



errath 2010; Bauer & Springel 2012), it has little effect on this par- 
ticular test. This owes in part to the implementation of the Balsara 
( 19891 switch in both cases, which reduces the artificial shear vis- 
cosity; without this, we see significantly greater damping of shear 
motions. 

In Fig. [8] we again repeat this experiment with the pressure- 
entropy formulation, but vary the smoothing kernel. We compare 
two spline kernels, for the neighbor numbers advocated in |Price| 
l |2012b| >; |Dehnen & Aly| ((2012), and the "core-triangle" kernels 
with sharp central peaks, proposed in Read et aL] {2010| >. For the 
standard initial conditions here, we see only subtle differences over 
a wide range from even the simplest kernel with 32 neighbors u 
through more than an order of magnitude higher neighbor number 
By explicitly eliminating the "surface tension" term, while main 
taining exact conservation, the pressure-entropy formulation ap 



9 Note that, at the same error level, the core-triangle kernels require signifi- 
cantly higher neighbor number, because they introduce a sharp bias towards 
the central particles. 



© 0000 RAS, MNRAS 000, 000-000 



10 Hopkins et al. 



Weak Instability 
Pressure- Entropy 




High Contrast 



Density-Entropy 



Pressure-Entropy 







ten*, ' J 








■Hi 







Figure 9. Comparison of Kelvin-Helmholtz behavior with different initial 
conditions in the Pressure-Entropy formulation at t 1 tkh (top) and f fS 
3tkh (bottom) on smoothing kernel (left: quintic spline with Nngb = 128; 
right: cubic spline with Wngb = 32). Compared to Fig. |6|8| we modify 
the initial conditions by multiplying the particle number, density, sound 
speed, pressure, shear velocity, and initial perturbation amplitude by fac- 
tors of = 0.5, 2.0, 2.0, 8.0, 0.5, 0.5, respectively. All of these changes in- 
crease the ratio of particle noise and pressure gradient errors relative to the 
KH growth. In this limit, the cubic spline with Wngb = 32, which involves 
larger gradient errors, is barely able to capture the instability; however the 
quintic spline with Nngb = 128 does well. "Standard" SPH completely fails 
to develop any "curls" in this limit. 



pears to significantly reduce the importance of kernel-level pressure 
gradient errors in obtaining the correct KH solution (compare e.g. 
the more significant kernel dependence found for standard but non- 
conservative density-entropy SPH in Read et al.|2010| Figs. 4 & 6). 
We have also tested the Wendland kernels proposed in |Dehnen &| 
|AlyH2012) and find (as they do) essentially identical performance 
to our spline results for the neighbor numbers here (although the 
kernels proposed there show greatly improved stability properties 
at higher neighbor number). 

However, in Fig. [9] we again examine the effects of different 
kernels (with the pressure-entropy formulation), but with different 
initial conditions. We reduce the particle number, shear velocity, 
and initial perturbation amplitude by factors of 2 and double the ini- 
tial sound speed; these changes all reduce the magnitude of the ini- 
tial KH instability and its growth rate, relative to the particle noise 
and error terms stemming from the kernel sum in the SPH pressure 
gradients that enter the EOM. This is designed to be challenging 
(even for grid codes). With this much weaker seeding and noisier 
particle distribution, we begin to see a dependence on the smooth- 
ing kernel. The simplest A^gb = 32 cubic spline, with the smallest 
neighbor number, leads to particle noise in the pressure gradients 
comparable to the actual signal. While the instability is still (barely) 
captured, the rolls are "too thin" and end up being sheared before 
they reach the proper height and can wrap appropriately. However, 
our standard quintic spline with A'ngb = 128 performs well, recov- 
ering all the key behaviors (note that the non-linear behavior is dif- 
ferent here than in the previous test, as it should be for the different 
pressure and shear velocities). Going to still higher resolution or 
varying the kernel at higher A'ngb gives well-converged results at 
this point. 

In Fig. [TO] we repeat our KH test again but multiply the initial 




Figure 10. Comparison of KH instabilities at t ?s 2 tkh (top) and t fti 1 tkh 
(bottom), for the "standard" (density-entropy) SPH (left) and pressure- 
entropy (right; with jc, = 1) SPH formulations. The initial conditions are 
as Fig. [6] but the initial density contrast is increased from a factor of 2 
to a factor of 20, with the initial particles no longer being equal mass but 
4 times more massive in the initial high-density region. Sharp boundary 
layers that lead to "gloopy" morphology are evident in standard SPH. In 
Pressure-Entropy SPH it remains well-behaved, although hints of "gloop- 
iness" appear in the transition between linear growth and fully non-linear 
instability. 



density contrast by a factor of 10, and instead of using constant- 
mass particles we use particles with masses a factor of 4 larger in 
the high-density region. As discussed in Read & Hayfield (20121, 
many proposed alternative formulations of SPH, designed to im- 
prove fluid mixing, fail in this regime (see e.g. Fig 6 in Re ad et al.| 
|2010]and Figs. Al & El in |Read & Hayfield|2012| >. And we see an 
even more pronounced "boundary layer" separating the phases in 
the standard density-entropy SPH formulation. This occurs because 
the higher density contrast exacerbates any (even small residual) 
surface tension term, and the multi-mass particles increase particle 
disorder and leading-order errors in the pressure gradient estimator. 
Multi-mass particles also make it critical to have a well-behaved 
criterion for smoothing lengths, and increase the errors from ne- 
glecting Vh terms. But we see that the pressure-entropy formu- 
lation remains well-behaved in this case. There is some increased 
hint of "gloopiness" as the instability transitions between linear and 
non-linear growth, but at least some of this is because of the (cor- 
rect) slower growth of the small-wavelength modes. 

Finally, it is worth noting (though perhaps as more of a cu- 
riosity) that the "poor" solution of the density-entropy formulation 
in the cases above looks very similar to the "correct" solution for 
a modestly magnetized medium (with some field component paral- 
lel to the shear or tangled; see e.g. prank et al.|1996| and references 
therein). This occurs because, if we consider the linear perturbation 
analysis of the KH instability, the "incorrect" surface tension force 
here is (for the initial linear stage) almost mathematically identical 
to a "correct" magnetic tension term for parallel fields with strength 

4.5 Rayleigh-Taylor Instabilities 

We now conside r the Raylei gh-Taylor (RT) instability, with initial 
conditions from 



Abel 



1 201 Ik. In a two-dimensional slice with 512 2 



particles and < x < 1/2 (periodic boundaries) and < y < 1 
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(reflecting boundary with particles at initial y < 0.1 or y > 0.9 
held fixed), we take j = 1.4 and initialize a density profile p(y) = 
Pi +(fl2 — pi)/(l+exp[— (y — 0.5)/A]) where pi = 1 and p 2 = 2 
are the density "below" and "above" the contact discontinuity and 
A = 0.025 is its width; initial entropies are assigned so the pres- 
sure gradient is in hydrostatic equilibrium with a uniform grav- 
itational acceleration g — — 1/2 in the y direction (at the inter- 
face, P = P2/7 = 10/7 so c s = 1). An initial y- velocity perturba- 
tion v y = 5v y (1 + cos (8-n- (x+ 1/4))) (1 + cos (5-7T (v -1/2))) with 
Sv y = 0.025 is applied in the range 0.3 < y < 0.7. 

Fig. [TT| shows the resulting evolution as a function of time, 
in the density-entropy and pressure-entropy formulations. As in 
SM12, both cases develop the RT instability with a similar linear 
growth time (slightly slower in the density-entropy case); we find 
that this is true for all of the kernel variations and both artificial 
viscosity choices discussed above, as well as for the slightly differ- 
ent (isoentropic) initial conditions used in SM12, also for 10-times 
smaller initial perturbations (5v y = 0.0025), and for resolutions as 
low as 50x100 particles[^However, they differ significantly in their 
nonlinear evolution; surface tension in the density-entropy formu- 
lation prevents the development of fine structure in the shear flow 
along the "fingers," obviously very closely related to the behavior 
in the KH tests. Pressure-entropy SPH, however, exhibits growth 
on the correct linear timescale and nonlinear behavior in agree- 
ment with that in Eulerian and moving-mesh schemes (e.g. Springel 
|20T0l >. 

4.6 The "Blob" Test 

We next consider the "blob" test, which is designed to test pro- 
cesses of astrophysical interest (e.g. ram-pressure stripping, mix- 
ing, KH and Rayleigh-Taylor instabilities in a multiphase medium). 
Again our initial conditions come from the Wengen test suite and 
are described in Age rtz et al.| ([2007). Briefly, the initial condi- 
tions include a spherical cloud of gas with uniform density in 
pressure equilibrium with an ambient medium, in a wind-tunnel 
with periodic boundary conditions. The imposed wind has Mach 
number M. = 2.7 and the initial density /temperature ratios are 
= 10. The box is a periodic rectangle with dimensions x, y, z = 
2000, 2000, 6000 kpc, with the cloud centered on 0, 0, -2000 kpc; 
the 9.6 x 10 6 particles are equal-mass and placed on a lattice. 

Fig.[T2]shows the resulting morphology of the cloud as a func- 
tion of time, in the pressure-entropy formulation. Fig. [T3] attempts 
to quantify the rate of cloud disruption following Agert z et al.| 
(20071, who defined a useful standard criterion for measuring the 
degree of mixing of initially dense blob gas: at each time, we sim- 
ply measure the total mass in gas with p > 0.64 p c and T < 0.9 T a 
(where p c and T a are the initial cloud density and ambient tempera- 
ture). We show the standard SPH (density-entropy) result from that 
paper, as compared to the prediction from the pressure-entropy for- 
mulation here. We also compare the result from the identical test in 
SMI 2, using the pressure-energy formulation but without includ- 
ing the V/z terms in the EOM, as well as the results from a high- 
resolution run with Enzo (O'Shea et al. 20041, a grid code. 

10 This agrees with the results in SM12, as well as the development of 
RT instabilities in blastwaves seen in other density-entropy implementa- 
tions (see e.g. Herant 1994 Price 2012c). It is not entirely clear why the 
RT instability fails to develop in standard SPH with the same initial con- 
ditions in Abel ( 201 1 1, but there are a number of additional differences in 
the algorithms (see §|5J. It may relate to the fact that the kernel instabilities 
discussed above can be much more severe in 2D standard SPH, requiring 
careful matching of neighbor number and kernel shape. 



The wind-cloud interaction generates a bow shock and imme- 
diately begins disrupting the cloud via a combination of KH and RT 
instabilities at its surface. By the middle panel in Fig.[T2] there is no 
visible large concentration of dense material. Compare this to the 
identical initial conditions run with standard SPH and grid codes in 
I Agertz et aL"1l2007| > (their Figs. 4 & 7). In standard SPH, the cloud 
is compressed to a "pancake," but the tension term prevents mixing 
at the surface and so a sizeable fraction survives disruption for large 
timescales. In contrast, the predicted morphology of the cloud here 
agrees very well with that in adaptive mesh codes therein (as well 
as moving-mesh methods in Sijacki et al. 2012). And quantitatively, 
we see this in Fig. 1 1 3 1 

However, once again we should note that the solution derived 
from the density-entropy formulation is remarkably similar to the 
"correct" solution with magnetic fields with /3 > 1 (compare e.g. 
|Mac Low et al.|1994]|Shin et aL]20 08 ), because the artificial surface 
tension term acts similarly to magnetic tension. 

4.7 Sub-Sonic Turbulence 

|Bauer & Sp ringel (2012) present a detailed study of the proper- 
ties of idealized driven isothermal turbulence in SPH - specifi- 
cally using GADGET-3 with the density-entropy formulation - as 
compared to moving mesh and grid methods (in the code AREPO; 
Springel 2010). Consistent with earlier results, they found that dif- 
ferent methods agree well when turbulence is super-sonic (e.g. Kit- 
|sionas et al.|2009| [Price & Federrath|2010) . But when the turbu- 
lence is sub-sonic, they found that density-entropy SPH tended to 
reproduce a smaller inertial range. However as discussed there and 
in Price (2012a), this can depend quite sensitively on the artificial 
viscosity prescriptions and other numerical details. Wc therefore re- 
produce the sub-sonic, driven turbulence experiment from Bauer & 
Springel (2012) but adopt the pressure-entropy formulation, to test 
whether it also differs significantly from the results in the density- 
entropy formulation. We adopt the identical setup and resolution 
(A. Bauer, private communication), with the initial conditions and 
driving algorithm described therein. Briefly, we initialize a box of 
unit length, density, and sound speed with 128 3 particles, and drive 
the turbulence in a narrow range of large-fc modes with character- 
istic Mach number M = 0.3 on the largest scales; unlike our other 
experiments, the gas is isothermal (7 = 1). We run the experiment 
to r = 25 (more than sufficient to reach steady-state). 

We compare the turbulent velocity spectrum measured follow- 
ing Bauer & Springel (2012) (we linearly interpolate the particle- 
centered velocity values to a uniform grid, fourier transform each 
velocity component, and average in bins of fixed k to obtain the 
kinetic energy power spectrum). The resulting power spectrum is 
plotted from k m 2n (though k < 10 are affected by the turbu- 
lent forcing as well as the finite box size) to k w 500 (a smooth- 
ing length h). We compare the pressure-entropy formulation to the 
results of standard (density-entropy formulation) SPH. The results 
are nearly identical (although there is slightly more power at in- 
termediate scales in the density-entropy formulation). The results 
follow the expected Kolmogorov inertial range, until ~ Ah, where 
they drop below owing to artificial numerical dissipation; they then 
rise at ~ 2h, owing to kemel-scale noise (largely from pressure 
gradient errors). 

We can compare this to other SPH results using the identi- 
cal initial conditions and driving routine. First, consider the re- 
sults from |Bauer & Springel] < |2012} , using "standard" GADGET-3 
(density-entropy EOM, without the additional algorithm improve- 
ments discussed in § [3} but w ith an artificial viscosity scheme fol- 
lowing Morris & Monaghan ( 1997) similar (though slightly differ- 
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Figure 11. Time evolution of a two-dimensional Rayleigh-Taylor (RT) instability (specific entropy map shown). Top: Standard (density-entropy) SPH. Bottom: 
Pressure-entropy formulation (.r, = 1). The RT instability develops in both cases but the mixing along the rising/sinking surfaces (a KH instability) is suppressed 
in density-entropy SPH, leading to different nonlinear outcomes. 



ent) to that adopted here, and the noisier Wngb = 32 cubic spline 
kernel. We also compare the results from [Pncel pOT2a) using a 
different SPH code (PHANTOM), with the same density-entropy 
EOM, and an artificial viscosity treatment identical to that here, 
but using an intermediate kernel (/Vngb = 58 cubic spline) and dif- 
ferent power spectrum calculation method. As shown there, the ar- 
tificial viscosity treatment dominates large/intermediate scales. So 
with identical artificial viscosity, the Price (2012a) result is identi- 
cal to ours over the inertial range (including the initial "drop off"). 
The smaller-scale (< 2 — 3 h) behavior is dominated by the kernel 
choice and method of power spectrum calculation. If we compare 
"standard" GADGET-3 with constant artificial viscosity I Gi ngold &| 
|Monaghan| 1983] > from [Bauer & Springel 1 201 2 ), we see much more 
severe excess dissipation, and almost no inertial range. 



Comparing to the results from AREPO in |Bauer & Spr ingel 
( |2012} (run either in fixed-grid or moving-mesh mode, they are 
nearly identical), we see that the deviations from Kolmogorov there 
are opposite in sign and quite distinct. In SPH, artificial viscos- 
ity produces excess dissipation on larger (resolved) scales, while 
kernel gradient errors lead to particle noise that boosts the power 
at the smallest scales. The former (larger-scale deviations), being 
dominated by artificial viscosity, are only indirectly tied to reso- 
lution; considering higher-resolution runs we confirm the results 
in Bauer & Springel (20121 regarding the relatively slow conver- 
gence the high-i spectrum in SPH relative to grid codes. The latter 
(smaller-scale deviations) are related to the conservative nature of 
the code and its ability to dissipate particle noise (their "return" to 
and overshoot of the Kolmogorov power-law are numerical, rather 
than physical consequences of a turbulent cascade). In Eulerian 
approaches, the lack of any but numerical viscosity concentrates 
the dissipation at the resolution/cell scale, which produces a "bot- 



tleneck" of excess power that cannot be dissipated at the larger 
(marginally resolved) scales; but this scales directly with the res- 
olution. 

4.8 Stellar Feedback in Isolated Galaxies 

Finally we test our modified algorithm on a fully non-linear sys- 
tem of particular astrophysical interest. Specifically, we re-run one 
of the star-forming galaxy disk models described in Hopkins et al. 
(|2011b| > and a subsequent series of papers ([Hopkins et al.|2012b|a| 
|201 la} . These simulations include gravity, collisionless stars and 
dark matter, and gas with a wide range of cooling processej^jfrom 
~ 10- 10 9 K, super-sonic turbulence, shocks, star formation in 
resolved giant molecular clouds, and explicit treatment of stellar 
feedback from SNe Types I & II, stellar winds, photoionization and 
photoelectric heating, and radiation pressure. Our purpose here is 
both to test how robust these algorithms are (since many problems 
with e.g. non-linear error terms or conservation errors often do not 
manifest in simple test problems such as those above), and to exam- 
ine how significant the pure numerical issues of fluid mixing are, 
relative to e.g. the previously-studied effects of adding/removing 
different physics in these simulations. 

We specifically re-run the "SMC" model (a model of an Small 
Magellanic Cloud-mass, gas-rich isolated dwarf galaxy), at "high" 
resolution (lpc softening length) as defined in the papers above, 
with all the above physical processes included exactly as in |Hop-| 
|kins et aT] ( |2012b"| >, but in one case adopting the density-entropy 
formulation and in another the pressure-entropy formulation. We 



" Cooling requires a density estimate, independent of whether one is 
needed for the EOM. For this, we use the standard SPH kernel density esti- 
mator, for the reasons in §13.41 
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Figure 13. Fraction of the initial cloud in Fig. [12] which remains both cold 
and dense (i.e. avoids mixing) as a function of time, relative to the KH 
growth time at the cloud surface. We compare the standard SPH (density- 
entropy) and pressure-entropy formulations here, as well as the pressure- 
energy formulation in SMI 2 (which does not include the V/i terms) and 
the results from a high-resolution grid code method (Enzo). The grid code 
and pressure formulations (independent of the V/i terms) agree reasonably 
well. Density formulations show slower mixing/stripping. 





Figure 12. Central density slices in the "blob test" (a high-density, pressure- 
equilibrium cloud hit by a wind) with the pressure-entropy formulation. 
Time (in units of the KH growth time) increases from top to bottom. The 
high-density (red/orange) cloud gas is efficiently mixed by instabilities 
within a couple cloud crossing times; the morphology and density distri- 
bution agree well with grid codes. 



choose this model because, of the galaxy types studied therein, it 
has the largest mass fraction in hot gas and is most strongly affected 
by thermal feedback mechanisms (e.g. SNe), and as a result is most 
sensitive to fluid mixing and phase structure. Figs. |15|T6| compare 
the visual morphologies and star formation rates from these sim- 
ulations. Taking into account that the systems are turbulent and 
non-linear, the morphologies are similar. They differ as one might 
expect: the pressure-entropy formulation increases mixing along 
phase boundaries, leading to less-sharp divisions between molec- 
ular clouds and/or hot, under-dense bubbles and the surrounding 
medium. This also makes the division between the star-forming 
and "extended" disk less sharp, as it is largely determined by the 
radius where the cooling rate becomes fast relative to the dynami- 
cal time[3 
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Figure 14. Velocity power spectrum in sub-sonic (M = 0.3) driven, 
isothermal turbulence. Each simulation uses 128 3 particles and identical 
driving. We compare the analytic Kolmogorov inertial-range model, to that 
calculated from our standard density-entropy (Eq. |14( and pressure-entropy 
(Eq. |21) formulations. The mean softening length h is shown as the verti- 
cal dotted line. Both do well down to a few softenings, where they first fall 
below the analytic result (excess dissipation) then rise above (kernel-scale 
noise). The EOM choice has a weak effect on the results. We compare dif- 
ferent SPH algorithms and codes (description in text); agreement is good 
where the same methods are used. The "PHANTOM" and "GADGET-3 
with Morris97 AV" use a density-entropy EOM with similar artificial vis- 
cosity treatment to our calculations (which dominates the inertial range), 
but different SPH kernels and power spectrum calculation methods (which 
dominate the noise at < 3h). The "standard" GADGET-3 calculation uses 
a constant artificial viscosity and different kernel, and produces almost no 
inertial range. The AREPO calculation considers a grid/moving mesh at the 
same resolution: deviations from Kolmogorov occur around the resolution 
limit, but are of a different character. 



12 We have also attempted a comparison using the algorithm in |Abel| 
(201 1) , from Fig. [2] Over short timescales this gives similar results to the 
pressure-entropy formulation. However, we cannot evolve the runs very 
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Figure 15. Morphology of the gas in a simulation of a star-forming disk 
galaxy (an isolated, SMC-mass dwarf") from Hopkins et al. ('2012b}, fol- 
lowing the resolved physics of the ISM and stellar feedback explicitly. 
Brightness encodes projected gas density (logarithmically scaled with 
6dex stretch); color encodes gas temperature with the blue material be- 
ing T < 1000K molecular and atomic gas, pink ~ 10 4 — 10 5 K warm ion- 
ized gas, and yellow > 10 6 K hot gas. We show the galaxy face-on top and 
edge-on bottom after several orbital periods of evolution (during which the 
galaxy is in quasi-steady state). We compare the identical simulation using 
the density-entropy (Eq. | 14) and pressure-entropy (Eq. |21) formulations of 
SPH. The two are largely similar. There is slightly sharper delineation be- 
tween phases (e.g. molecular clouds and hot bubbles) in the density-entropy 
formulation; this includes the sharper transition between the star-forming 
and outer disks (physically, where the cooling time becomes longer than 
the dynamical time). 

Quantitatively, the time-averaged star formation rates differ 
remarkably little, just w 20% (with higher SFRs in the pressure- 
entropy case); this owes to the fact that they are globally set by 
equilibrium between momentum injection from stellar feedback 
and dissipation of turbulent energy/momentum on a dynamical time 
(Hopk ins et al.|20 12b). The instantaneous rates vary more rapidly 
and can differ by factors of ~ 2 — 3. We also compare the mass- 
loading of galactic winds in each simulation, defined as in Hopkins 
|et al.| ( [2012a) as the ratio of total "wind mass" (material with posi- 
tive Bernoulli parameter, i.e. which will be unbound in the absence 
of external pressure forces) to new stellar mass formed since the 
beginning of the simulation. Here the qualitative behavior is the 
same but the wind mass-loading is systematically smaller in the 
pressure-entropy case, by a factor of ~ 1.5 — 2. This is directly re- 
lated to the higher SFR; both owe to the increased mixing reducing 
the cooling time of the hot gas bubbles (which provide, for this 
small galaxy, much of the wind mass). In more massive galaxies, 
the winds in these simulations are less dominated by hot gas, so the 
effect will be smaller. These differences (while not negligible) are 
smaller than those typically caused by adding or removing differ- 
long before a situation similar to Fig. [2] arises when, for example, a SNe 
occurs in an under-dense region and its energy is deposited in a small num- 
ber of particles, leading to large conservation errors. 
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Figure 16. Star formation rate and stellar wind mass-loading (ratio of total 
unbound gas mass to total new stellar mass formed since the start of the 
simulation) as a function of time, for the galaxy simulations in Fig. |15| 
Although stochastic variations in the SFR differ by factors ~ 2 — 3, the 
time-averaged SFR is only ~ 20% lower in the density-entropy formulation. 
The wind mass-loading is systematically higher in this case by a factor ~ 
1.5 — 2.0. This stems from increased phase mixing in the pressure-entropy 
formulation introducing more cold, dense gas into "hot" gas, enhancing its 
cooling before it can vent out of the disk. 

ent feedback mechanisms, or even changing details of their imple- 
mentation (shown in Hopkins et al. 2012a). For example, removing 
stellar feedback entirely in this simulation changes the SFR by a 
factor of ~ 50! Integrated over cosmological times, or allowing for 
different galaxy conditions because of the integrated effects of dif- 
ferent heating/cooling efficiencies as gas actually accretes into ha- 
los, however, these small differences can easily build up into more 
significant divergence (see e.g.JVogelsberger et al .|201 1} . 

5 DISCUSSION 

With inspiration from SMI 2, we re-derive a fully self-consistent set 
of SPH equations of motion which are independent of the kernel- 
calculated density, and therefore remove the well-known "surface 
tension" terms that can suppress fluid mixing. The equations still 
depend on the medium having a differentiable pressure (hence still 
require some artificial viscosity to capture shocks), but unlike the 
traditional SPH EOM, remain valid through contact discontinuities. 

Our derivation of the EOM relies on the key conceptual point 
in SMI 2, that the SPH volume element does not have to explic- 
itly involve the mass density. However, our derivation and resulting 
EOM adds four key improvements: (1) We rigorously derive the 
equations from the discretized particle Lagrangian. This guaran- 
tees one of the most powerful features of SPH, namely manifest 
simultaneous conservation of energy, entropy, momentum, and an- 
gular momentum, and an exact solution to the particle continuity 
equation. (2) We similarly derive the "V/t" terms, which are re- 
quired for manifest conservation if the SPH smoothing lengths hi 
are not everywhere constant. (3) We derive an "entropy formula- 
tion" of the equations that allows for the direct evolution of the 
entropy, avoiding the need to construct/evolve an energy equation, 
and gives better entropy conservation properties as in S02; this also 
happens to minimize the correction terms involved in using a "par- 
ticle neighbor number" definition to define h, as compared to the 
"energy formulation." (4) We show how the Lagrangian derivation 
can be generalized to separate definitions of the thermodynamic 
volume element (relating e.g. P and u) and that used to define the 
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smoothing lengths. This resolves problems of numerical stability 
and excess diffusion in strong shocks and/or large density contrasts, 
and automatically allows for varying particle masses. 

In fact, we derive a completely general, Lagrangian form of 
the EOM, including the Vh terms, for any definition of the SPH 
thermodynamic volume element. Essentially any particle-carried 
quantity can be used in the kernel sum entering the EOM, and 
any (not necessarily the same) differentiable function used to define 
how the smoothing lengths hi are scaled. In some ways this replaces 
the long-known "free weighting functions" used to define the SPH 
EOM in their original "discretized volume element" formulation. 
However, in that approach, the choice of different functions gener- 
ically violates conservation and continuity; here, we demonstrate 
that a similar physical degree of freedom can be utilized in the dis- 
cretization of the equations of motion without such violations. 

Based on this degree of freedom, it is easy to see how different 
discretizations of the EOM might be optimized for some problems. 
By choosing the required kernel-evaluated element to directly rep- 
resent a very smooth/stable property in the system, one not only re- 
moves spurious "tension" terms associated with discontinuities in 
other system variables, but also minimizes the inevitable discretiza- 
tion error from representing these quantities with a kernel sum. 
For the constant-pressure (but mixed density) fluid mixing tests 
we show here, the optimal choice is the "pressure formulation." 
In MHD applications this kernel sum could trivially be altered to 
include the magnetic pressure. However if simulating an incom- 
pressible or weakly compressible fluid, the "density formulation" 
may well be superior. Direct kernel sums of nominally "higher or- 
der" properties such as the vorticity or vortensity are also valid and 
may represent useful formulations for some problems. It is even 
possible (in principle) to generalize our derivation to one in which 
different particle subsets have differently-defined volume elements; 
although we caution that such an approach requires great care. 

For the test problems here, we show that the "pressure- 
entropy" and "pressure-energy" formulations dramatically improve 
the treatment of fluid interface instabilities including the Kelvin- 
Helmholtz instability, Rayleigh-Taylor instability, and the "blob 
test" (a mix of Kelvin-Helmholtz and Raleigh-Taylor instabilities 
as well including non-linear evolution and shock capturing); giving 
results very similar to grid methods. They also remove the "deform- 
ing" effect of the surface tension term (allowing, for example, the 
long-term evolution of an irregular shape of gas at constant pressure 
but high density contrast); deformation is difficult to avoid even 
in grid codes (unless the chosen geometry matches the grid), and 
would otherwise require moving-mesh approaches to follow. How- 
ever, unlike some of the modifications in the literature proposed 
to improve the fluid mixing in SPH (which violate conservation), 
the manifest conservation properties of our derivation mean that it 
remains well-behaved even in very strong shocks and does not en- 
counter problems of either energy conservation or particle order in 
e.g. extremely strong blastwave problems. 

With these changes in place, we find weaker (albeit still signif- 
icant) residual effects from improvement in the artificial viscosity 
scheme. Comparisons of such schemes are well-studied and what 
we implement here is still not the most sophisticated possible treat- 
ment, although it still considerably reduces artificial viscosity away 
from shocks (for more detailed studies, see e.g. |Cullen & Dehnen| 
|2010[ >. We find similar effects from changes to the SPH smoothing 
kernel. Our favored kernel is taken from more detailed kernel com- 
parison studies in Hon gbin & X?ii| ( [2005) 1; [Dehnen & Aly] pOT2] l; 
however, unlike some other SPH formulations, we find that even 
the "simplest" kernel possible (/Vngb = 32 cubic spline) reproduces 



good results in several tests, except in the expected regime where 
we wish to resolve kernel-scale growing instabilities that rely on 
sub-sonic motions at the level of M < A'ngb an d so summation er- 
rors dominate. This relates to the manifest conservation and mainte- 
nance of good particle order implicit in the EOM (see Price 2012b| >. 

We test the algorithm not just in the "standard" set of test prob- 
lems but also an example of direct astrophysical interest, simulat- 
ing the evolution of galaxies with a multi-phase ISM. This is useful 
because it makes clear that for this problem, at least, the differ- 
ences arising from the treatment of different physics (e.g. how cool- 
ing, star formation, stellar feedback, and AGN feedback are imple- 
mented) makes, on average, larger differences than the numerical 
scheme (also shown in other code comparisons; e.g. |Scannapieco| 
|et al.| [2012). This is not surprising: those choices lead to orders- 
of-magnitude differences as opposed to the (still significant) factor 
~couple effects of numerical choices. Moreover the differences we 
are concerned with here largely pertain to mixing in sub-sonic, non- 
radiative flows dominated by thermal pressure; in contrast many as- 
trophysical problems of interest involve highly super-sonic, radia- 
tive, gravity-dominated flows. In that limit, the differences owing 
to the algorithm are often - though certainly not always - mini- 
mized (see references in §[TJ. But there are important regimes with 
transonic flows where the numerical approach can make larger dif- 
ferences (e.g. cosmological inflows & outflows; see Vogelsberger 
|et al.|201 l|[Keres et al.|2"0l2"{|Torrey et al.|201 1| >. Even in idealized 
test problems, we caution that simple physical differences can pro- 
duce larger distinctions than the numerical method. For example, 
for several fluid mixing problems considered here, the "correct" 
MHD solution in the presence of an equipartition magnetic field 
can resemble the "standard" (density-entropy) SPH solution with- 
out a magnetic field, as opposed to the results from our pressure- 
entropy formulation or grid codes without such fields. The reason 
is that the real magnetic tension suppresses mixing, similar to the 
(purely numerical) "surface tension" term discussed in the text (so 
one might obtain a more "realistic" solution, but for entirely wrong 
reasons). 

Ultimately, the numerical formulations derived here should 
provide the basis for a more rigorous approach to the "flavors" of 
SPH, and a means to compare the consequences of the fundamental 
choice of how to discretize any SPH approach. This change to the 
algorithm is not a panacea! Fortunately, the modified equations of 
motion proposed here can be trivially incorporated with many other 
methods that improve on other numerical aspects, for example the 
inviscid algorithm in Cullen & Deh nen] j2010) , the higher-order 
dissipation switches in |Price| j2008) ; |Rosswog| ( |2010^ and [Read] 
|& Hay field (2012), and/or the gradient error reducing integral for- 
mulation of the kernel equations in Garcfa-Senz et al. (2012). We 
wish to stress that - although SPH certainly has some disadvan- 
tages which we have not attempted to address here - poor fluid 
mixing in contact discontinuities is not necessarily an "inherent" 
property of SPH. This problem can be improved without requiring 
additional dissipation terms (and without additional computational 
expense) while retaining what is probably the greatest advantage of 
SPH algorithms, namely their excellent conservation properties. 
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